#Figure A1: The Causal Effects of Policy Features on Long-Term Investment by Subgroups
#install.packages("foreign")
#install.packages("ggplot2")
#install.packages("dplyr")
#install.packages("xlsx")
#install.packages("gridExtra")

library(foreign)
library(ggplot2)
library(dplyr)
library(xlsx)
library(gridExtra)

#Set working directory
#setwd("C:/Users/mmann/Dropbox/Disaster policy/02 Empirics/10 Replicationarchive/03 Results")
d = read.xlsx("figA1_res_diff.xlsx", sheetIndex = 1)
hl = read.xlsx("figA1_res_est.xlsx", sheetIndex = 1)

#Prepare data
#---------------------------
d = rename(d, trust = NA., trust_pe = y1, trust_se = r1,
           respdis = NA..1, respdis_pe = y1.1, respdis_se = r1.1,
           income = NA..2, income_pe = y1.2, income_se = r1.2,
           sex = NA..3, sex_pe = y1.3, sex_se = r1.3,
           uni = NA..4, uni_pe = y1.4, uni_se = r1.4,
           age = NA..5, age_pe = y1.5, age_se = r1.5)

d = d[c(1:3, 9, 11), ] #extract only data needed

hl = filter(hl, NA. == "2.framing" | NA. == "3.framing")
hl = select(hl, 19:36) 
hl = rename(hl, trust = NA..6, trust_pe = y1.6, trust_se = r1.6,
           respdis = NA..7, respdis_pe = y1.7, respdis_se = r1.7,
           income = NA..8, income_pe = y1.8, income_se = r1.8,
           sex = NA..9, sex_pe = y1.9, sex_se = r1.9,
           uni = NA..10, uni_pe = y1.10, uni_se = r1.10,
           age = NA..11, age_pe = y1.11, age_se = r1.11)

#Political Trust
trust = select(hl, 1:3)
trust$grouptrust = "Trust: High"
trust = select(trust, -trust)
newrow = c(0,0,"Trust: High")
trust = rbind(trust,newrow)

trust1 = select(d, 1:3)
trust1$grouptrust = ifelse(trust1$trust=="2.framing#1.trust_high" | trust1$trust=="3.framing#1.trust_high","Difference","Trust: Low")
trust1 = select(trust1, -trust)

trust = rbind(trust,trust1)
var = as.data.frame(c("Prep", "Prep", "Baseline", "Baseline", "relief", "relief", "diff_prep", "diff_relief"))
names(var) = "var"
trust = cbind(var,trust)
trust$trust_pe = as.numeric(trust$trust_pe)
trust$trust_se = as.numeric(trust$trust_se)

#Responsibility Beliefs
respdis = select(hl, 4:6)
respdis$grouprespdis = "Prep. Responsbility: Gvt."
respdis = select(respdis, -respdis)
newrow = c(0,0,"Prep. Responsbility: Gvt.")
respdis = rbind(respdis,newrow)

respdis1 = select(d, 4:6)
respdis1$grouprespdis = ifelse(respdis1$respdis=="2.framing#1.respdis_high" | respdis1$respdis=="3.framing#1.respdis_high","Difference","Prep. Responsibility: Each Person")
respdis1 = select(respdis1, -respdis)

respdis = rbind(respdis,respdis1)
respdis$respdis_pe = as.numeric(respdis$respdis_pe)
respdis$respdis_se = as.numeric(respdis$respdis_se)

#Income
income = select(hl, 7:9)
income$groupincome = "Income: High"
income = select(income, -income)
newrow = c(0,0,"Income: High")
income = rbind(income,newrow)

income1 = select(d, 7:9)
income1$groupincome = ifelse(income1$income=="2.framing#1.income_high" | income1$income=="3.framing#1.income_high","Difference","Income: Low")
income1 = select(income1, -income)

income = rbind(income,income1)
income$income_pe = as.numeric(income$income_pe)
income$income_se = as.numeric(income$income_se)

#Gender
sex = select(hl, 10:12)
sex$groupsex = "Sex: Male"
sex = select(sex, -sex)
newrow = c(0,0,"Sex: Male")
sex = rbind(sex,newrow)

sex1 = select(d, 10:12)
sex1$groupsex = ifelse(sex1$sex=="2.framing#1.sex" | sex1$sex=="3.framing#1.sex","Difference","Sex: Female")
sex1 = select(sex1, -sex)

sex = rbind(sex,sex1)
sex$sex_pe = as.numeric(sex$sex_pe)
sex$sex_se = as.numeric(sex$sex_se)

ordering_gender = data.frame(c(5,6,4,3,1,2,7,8))
names(ordering_gender) = "ordering"
sex = cbind(sex, ordering_gender)
sex = arrange(sex, sex$ordering)
sex = select(sex, -ordering)

#Education
uni = select(hl, 13:15)
uni$groupuni = "University Degree: Yes"
uni = select(uni, -uni)
newrow = c(0,0,"University Degree: Yes")
uni = rbind(uni,newrow)

uni1 = select(d, 13:15)
uni1$groupuni = ifelse(uni1$uni=="2.framing#1.uni_degree" | uni1$uni=="3.framing#1.uni_degree","Difference","University Degree: No")
uni1 = select(uni1, -uni)

uni = rbind(uni,uni1)
uni$uni_pe = as.numeric(uni$uni_pe)
uni$uni_se = as.numeric(uni$uni_se)

#Age
age = select(hl, 16:18)
age$groupage = "Age: Old"
age = select(age, -age)
newrow = c(0,0,"Age: Old")
age = rbind(age,newrow)

age1 = select(d, 16:18)
age1$groupage = ifelse(age1$age=="2.framing#1.age_high" | age1$age=="3.framing#1.age_high","Difference","Age: Young")
age1 = select(age1, -age)

age = rbind(age1, age)
age$age_pe = as.numeric(age$age_pe)
age$age_se = as.numeric(age$age_se)

ordering_age = data.frame(c(4,5,6,7,8,1,2,3))
names(ordering_age) = "ordering"
age = cbind(age,ordering_age)
age = arrange(age, age$ordering)
age = select(age, -ordering)

#Put all of the above together and transform data
d = cbind(trust,respdis,income,sex,uni,age)

sqrf <- function(x){
  sqrt(x)
}

dsub <- select(d, trust_se, respdis_se, income_se, sex_se, uni_se, age_se)


dse <- lapply(dsub, sqrf)
dse <- as.data.frame(dse)
d = select(d, -trust_se, -respdis_se, -income_se, -sex_se, -uni_se, -age_se)
d <- cbind(d, dse)

ordering = data.frame(c(4,6,2,1,3,5,7,8))
names(ordering) = "ordering"
d = cbind(d,ordering)
d = arrange(d, d$ordering)


#Generating upper and lower bounds and preparing figure
#Political Trust
d$trust_lower <- d$trust_pe  - 1.96*d$trust_se
d$trust_upper <- d$trust_pe  + 1.96*d$trust_se
order = data.frame(c(0.8,1.2,1.8,2.2,3.8,4.2,3,5))
names(order) = "order"
d = cbind(d,order)
d$order <- d$order * -1 + 6

d$grouptrust <- as.factor(d$grouptrust)
d$grouptrust <- factor(d$grouptrust, levels = c(levels(d$grouptrust)[3],levels(d$grouptrust)[2],levels(d$grouptrust)[1]))

#Responsibility Beliefs
d$respdis_lower <- d$respdis_pe  - 1.96*d$respdis_se
d$respdis_upper <- d$respdis_pe  + 1.96*d$respdis_se

d$grouprespdis <- as.factor(d$grouprespdis)
d$grouprespdis <- factor(d$grouprespdis, levels = c(levels(d$grouprespdis)[3],levels(d$grouprespdis)[2],levels(d$grouprespdis)[1]))

#Income
d$income_lower <- d$income_pe  - 1.96*d$income_se
d$income_upper <- d$income_pe  + 1.96*d$income_se

d$groupincome <- as.factor(d$groupincome)
d$groupincome <- factor(d$groupincome, levels = c(levels(d$groupincome)[3],levels(d$groupincome)[2],levels(d$groupincome)[1]))

#Gender
d$sex_lower <- d$sex_pe  - 1.96*d$sex_se
d$sex_upper <- d$sex_pe  + 1.96*d$sex_se

d$groupsex <- as.factor(d$groupsex)
d$groupsex <- factor(d$groupsex, levels = c(levels(d$groupsex)[2],levels(d$groupsex)[3],levels(d$groupsex)[1]))

#Education
d$uni_lower <- d$uni_pe  - 1.96*d$uni_se
d$uni_upper <- d$uni_pe  + 1.96*d$uni_se

d$groupuni <- as.factor(d$groupuni)
d$groupuni <- factor(d$groupuni, levels = c(levels(d$groupuni)[2],levels(d$groupuni)[3],levels(d$groupuni)[1]))

#Age
d$age_lower <- d$age_pe  - 1.96*d$age_se
d$age_upper <- d$age_pe  + 1.96*d$age_se

d$groupage <- as.factor(d$groupage)
d$groupage <- factor(d$groupage, levels = c(levels(d$groupage)[2],levels(d$groupage)[1],levels(d$groupage)[3]))

#Add labels
labs <- c(expression(paste(Delta, " Compensation")), "Compensation", expression(paste(Delta, " Damage Reduction")), 
          "Damage Reduction", "Control") 

#Create figures
figa1_trust <- ggplot(d, aes(x=order, y=trust_pe, colour=grouptrust)) + 
  geom_point(aes(x = order, y = trust_pe, colour = grouptrust), show.legend = T, size = 2) +
  geom_segment(aes(x = order, xend = order, y = trust_lower, yend = trust_upper), show.legend = T, size=1) +
  coord_flip(ylim = c(-8, 10)) + theme_bw() +
  scale_x_continuous(name="", breaks=seq(1,5,1),labels=labs) +
  scale_y_continuous(name="Change in Long-Term Investment (million $)",
                     breaks = c(-8,-6,-4,-2,0,2,4,6,8,10)) +
  geom_hline(yintercept = 0,size=0.5,colour="black",linetype="dashed") +
  theme(legend.position = c(.85, 0.91),
        legend.background = element_rect(fill="transparent"),
        legend.key = element_blank(), axis.text = element_text(size = 16), axis.title.x=element_text(size=18),
        legend.text=element_text(size=13)) +
  scale_colour_grey(name = "") 
ggsave("FigureA1a_trust.pdf", width=10, height=5)


figa1_respbeliefs <- ggplot(d, aes(x=order, y=respdis_pe, colour=grouprespdis)) + 
  geom_point(aes(x = order, y = respdis_pe, colour = grouprespdis), show.legend = T, size = 2) +
  geom_segment(aes(x = order, xend = order, y = respdis_lower, yend = respdis_upper), show.legend = T, size=1) +
  coord_flip(ylim = c(-6, 10)) + theme_bw() +
  scale_x_continuous(name="", breaks=seq(1,5,1),labels=labs) +
  scale_y_continuous(name="Change in Long-Term Investment (million $)",
                     breaks = c(-6,-4,-2,0,2,4,6,8,10)) +
  geom_hline(yintercept = 0,size=.5,colour="black",linetype="dashed") +
  theme(legend.position = c(.85, 0.91),
        legend.background = element_rect(fill="transparent"),
        legend.key = element_blank(), axis.text = element_text(size = 16), axis.title = element_text(size = 16)) +
  scale_colour_grey(name = "") 
ggsave("FigureA1b_respbeliefs.pdf", width=10, height=5)

fig2_income <- ggplot(d, aes(x=order, y=income_pe, colour=groupincome)) + 
  geom_point(aes(x = order, y = income_pe, colour = groupincome), show.legend = T, size = 2) +
  geom_segment(aes(x = order, xend = order, y = income_lower, yend = income_upper), show.legend = T, size=1) +
  coord_flip(ylim = c(-6, 10)) + theme_bw() +
  scale_x_continuous(name="", breaks=seq(1,5,1),labels=labs) +
  scale_y_continuous(name="Change in Long-Term Investment (million $)",
                     breaks = c(-6,-4,-2,0,2,4,6,8,10)) +
  geom_hline(yintercept = 0,size=.5,colour="black",linetype="dashed") +
  theme(legend.position = c(.85, 0.91),
        legend.background = element_rect(fill="transparent"),
        legend.key = element_blank(), axis.text = element_text(size = 16), axis.title = element_text(size = 16)) +
  scale_colour_grey(name = "") 
ggsave("FigureA1c_income.pdf", width=10, height=5)

figa1_gender <- ggplot(d, aes(x=order, y=sex_pe, colour=groupsex)) + 
  geom_point(aes(x = order, y = sex_pe, colour = groupsex), show.legend = T, size = 2) +
  geom_segment(aes(x = order, xend = order, y = sex_lower, yend = sex_upper), show.legend = T, size=1) +
  coord_flip(ylim = c(-6, 10)) + theme_bw() +
  scale_x_continuous(name="", breaks=seq(1,5,1),labels=labs) +
  scale_y_continuous(name="Change in Long-Term Investment (million $)",
                     breaks = c(-6,-4,-2,0,2,4,6,8,10)) +
  geom_hline(yintercept = 0,size=.5,colour="black",linetype="dashed") +
  theme(legend.position = c(.85, 0.91),
        legend.background = element_rect(fill="transparent"),
        legend.key = element_blank(), axis.text = element_text(size = 16), axis.title = element_text(size = 16)) +
  scale_colour_grey(name = "") 
ggsave("FigureA1d_gender.pdf", width=10, height=5)

figa1_educ <- ggplot(d, aes(x=order, y=uni_pe, colour=groupuni)) + 
  geom_point(aes(x = order, y = uni_pe, colour = groupuni), show.legend = T, size = 2) +
  geom_segment(aes(x = order, xend = order, y = uni_lower, yend = uni_upper), show.legend = T, size=1) +
  coord_flip(ylim = c(-6, 10)) + theme_bw() +
  scale_x_continuous(name="", breaks=seq(1,5,1),labels=labs) +
  scale_y_continuous(name="Change in Long-Term Investment (million $)",
                     breaks = c(-6,-4,-2,0,2,4,6,8,10)) +
  geom_hline(yintercept = 0,size=.5,colour="black",linetype="dashed") +
  theme(legend.position = c(.85, 0.91),
        legend.background = element_rect(fill="transparent"),
        legend.key = element_blank(), axis.text = element_text(size = 16), axis.title = element_text(size = 16)) +
  scale_colour_grey(name = "") 
ggsave("FigureA1e_educ.pdf", width=10, height=5)

fig1_age <- ggplot(d, aes(x=order, y=age_pe, colour=groupage)) + 
  geom_point(aes(x = order, y = age_pe, colour = groupage), show.legend = T, size = 2) +
  geom_segment(aes(x = order, xend = order, y = age_lower, yend = age_upper), show.legend = T, size=1) +
  coord_flip(ylim = c(-8, 10)) + theme_bw() +
  scale_x_continuous(name="", breaks=seq(1,5,1),labels=labs) +
  scale_y_continuous(name="Change in Long-Term Investment (million $)",
                     breaks = c(-8,-6,-4,-2,0,2,4,6,8,10)) +
  geom_hline(yintercept = 0,size=.5,colour="black",linetype="dashed") +
  theme(legend.position = c(.85, 0.91),
        legend.background = element_rect(fill="transparent"),
        legend.key = element_blank(), axis.text = element_text(size = 16), axis.title = element_text(size = 16)) +
  scale_colour_grey(name = "") 
ggsave("FigureA1f_age.pdf", width=10, height=5)







